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Abstract 

We present an algorithm for computing the main topofogical charac- 
teristics of three-dimensional bodies. The algorithm is based on a dis- 
cretization of Morse theory and uses discrete analogs of smooth func- 
tions with only nondegenerate (Morse) and the simplest degenerate 
critical points. 
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1 Introduction 

The problem of computing topological characteristics of geometrical objects 
of complicated shape appears in many applications. This led to the rise 
of computational topology which is fast developing now [H El [3]. Usually 
a problem of computational topology reduces to effective computation of 
Betti numbers, i.e. the ranks of homology groups. 

One of the most typical approaches to computing homology groups is 
based on Morse theory [H [5] which has to be modified for studying objects 
that are defined by discrete data. Discrete approaches to the Morse theory 
were systematically studied in [T] and [6] and various discrete realizations of 
Morse theory were considered in [71 El O [10] . 

In this article we propose a method for computing the Betti numbers of 
three-dimensional bodies. Therewith the bodies are defined as unions of unit 
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cubes with integer-valued coordinate of vertices in the three-dimensional 
Euclidean space. The main feature of the approach is 

• the use of discrete analogs of smooth functions (of two variables) which 
in addition to standard nondegenerate (Morse) critical points have the 
simplest degenerate critical points of the type of the '^monkey saddle" 
in the case under study. 

The scheme of this approach is as follows. Let M be a three-dimensional 
body with a boundary in the three-dimensional Euclidean space. We con- 
sider the critical points of a differentiable (smooth) function / on M, i.e. the 
points at which the gradient of / vanishes. The values of / at these points 
are called critical. When the value of / passes through a critical value the 
topology of the excursion set {/ < const } changes. By using critical points 
we construct a chain complex C* = Cq © Ci © C2 which is a graded com- 
mutative group. To the critical points there correspond the basic elements 
(generators) of the complex C\ as of a vector space over the commutative 
field Z2 which consists of two elements. To a nondegenerate critical point 
of index i there corresponds a generator of Cj and to a critical point of the 
"monkey saddle" type there correspond two generators of Ci. The gradi- 
ent flow of / defines the boundary operators which arc homomorphisms of 
commutative groups di : Ci ^ Cj-i. The final part of this algorithm is 
standard: since the kernel Ker5 of a homomorphism 5 : ^ C* lies in 
its image Im5 the homology groups (with coefficients in Z2) are correctly 
defined as the quotient spaces 

^ ^ Kcr(9, ^Q-i) 
lni{di+i : Ci+i ^ Ci)' 

The dimension of the i-th homology group as of a vector space over Z2 
is called the i-th. Betti number of M and it is denoted as 6,. The Betti 
numbers bo, bi, and b2 which are interpreted as the number of connected 
components, the number of linearly independent one-dimensional cycles and 
the number of holes in M and their alternated sum x = ^0 — ^1 + ^2 are the 
main topological invariants of M. 

Clearly, we consider a discrete version of this algorithm. Note that since a 
body is three-dimensional it is enough to restrict consideration to homology 
with Z2 coefficients for a complete description of homological characteristics 
of M. Also we may use the Alexander duality and compute only the homo- 
morphism 9i : Ci — )• Cq. We will discuss that here in morein detail. 

Since we consider only bounded bodies, every such a body is assumed to 
be embedded into a cube K in the three-dimensional space. The algorithm 
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of computing the matrices of boundary operators has complexity 0{n) where 
n is the number of vertices, points with integer-valued coordinates, in K, 
and the complexity of computing Betti numbers is equal to O(n^), where ric 
is the number of critical points of /. 

Remark 1. There exist infinitely many different topological types of 
critical points of analytic functions of n (n > 2) variables and their classifi- 
cation is very complicated. However there are finitely many Morse singular- 
ities and topological surgeries of corresponding level surfaces are completely 
well described. Therefore in computational topology it was supposed to 
use only the functions with Morse singularities for computing topological 
characteristics (see the scheme below). However it appears that after dis- 
cretization the classification of critical points becomes finite which is quite 
natural and by describing all topological surgeries corresponding to non- 
Morse singularities one may use more general functions for realization of the 
scheme mentioned above and we demonstrate that in the three-dimensional 
case. 

Remark 2. There are algorithms for computing Betti numbers with 
complexity asymptotically linear in n,[l] which substantially use the specifics 
of the three-dimensional case (see, for instance, [11)). In applications [12] 
in which our algorithm was used the number ric of critical points is small 
as compared with n and that gives a practical advantage. Moreover the 
proposed method may be generalized to higher dimensions. 

2 Critical points of the diagonal function on a Eu- 
clidean cube 

By an elementary interval / C M we mean the set / = [/,/ + 1] for some 
/ € Z. We analogously define the elementary square Q = /i x /2 C 
and the elementary cube C = Ii x I2 x c , where 1^ are elementary 
intervals. 

We consider three-dimensional bodies M formed by elementary cubes 
and lying inside the bounded domain 

K = [0, N] X [0, N] X [0, N] C M^ 

for some natural number A^. 

^The complexity of efTective algorithms of distinguishing connected components is equal 
to 0{na{n)) where a{n) is a slowly growing function which is inverse to the the Ackermann 
function. 
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In fact we compute homology not of M but of the body M which is 
constructed from M by unstacking joint vertices and boundary edges for 
such elementary cubes, in M, which intersect only at these joint vertices 
and edges. 

We make to remarks: 

1) we consider M due to applications of the algorithm; H 

2) the body M itself, or more precisely a topologically equivalent body, 
may be realized by unions of elementary cubes in a cube K' of a larger size, 
as shown by the following construction. Let us consider the cube 

K' = [0,3N] X [0,3iV] X [0,3iV], 

and the linear map P{xi, X2, X3) = (3xi, 8x2, 8x3) from K to K' . We put 
M' = P{M), i.e.. M' is a subdivision of M such that every elementary 
edges splits into three parts and every elementary cube splits into 9 parts. 
Then we remove from m' all elementary cubes such that they have a joint 
face with dM' . We obtain the new cubical set M. In the sequel we assume 
that M is obtained, if need be by such a surgery. 

Let us take for a function f on K the "diagonal" function 

/(Xl, X2, X3) = Xl + X2 + X3. 

As usual, we put 

M"" = {x G M|/(x) < a}. 

It is clear that if oi < 02 and the set /~^([ai,a2]) fl M does not contain 
points with integer-valued coordinates, then M^^ and M"^ are topologically 
equivalent (see, for instance, [3]). Therefore, under a continuous variation 
of a from to 3A^ the topological type of M"" changes only when a passes 
through a value of the form a = /(xi, X2, X3) where Xj € Z. 

Let ki, k2, be a triple of integer numbers such that v = {ki,k2, k^) € 
K. We put 

N{v) = {x£ M\\xi -ki\<l,i = 1,2,3}. 

Hence N{v) consists of eight elementary cubes surrounding v and lying in 
M. For a = f{v) and |6 — a| < g we put 

NyM'' = N{v) n M\ 

In [12] this algorithm is used for evaluation of topological characteristics of oil reser- 
voirs. Such a reservoir is modeled by a three-dimensional body that is a union of three- 
dimensional cells. Therewith it is assumed that a passing of oil from one cell to another 
is possible only though the joint two-dimensional face, and there is no passing through a 
joint vertex or a joint edge. To this end, we have to make a preliminary processing of M. 
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We say that v £ M is a. critical point of a function / if for sufficiently 
small positive e < | the sets NyM"'~'^ and A^'^M""'"^ are not topologically 
equivalent. 

For distinguishing the critical point type we code the cubes from N{v) 
as follows. This neighborhood consists of eight elementary cubes and with 
respect to the central integer-valued point every cube may be defined by 
changes of the coordinates whose values are either greater or less inside the 
cube than at v. Hence each cube is coded by a triple of signs "± it ±". 
Every such a cube may lie or not lie in M. Hence we may correspond to 
N{v) the pair of matrices 



where t*** = 1 if the cube lies in M and t**,, = otherwise. 

We say that a critical point w G M is nondegenerate if NyM°'~^^ topologi- 
cally looks like the body N^M'^~^ with a handle of index z, i = 0, 1, 2, glued 
to its boundary, i.e. one of the three possibilities holds: 

1) (a handle of index 0) NyM"'~^ is empty and there appears a new 
connected component N^M""^^; 

2) (a handle of index 1) neighborhoods of two points on the boundary 
of NyM'^^^ are connected by a thickened interval; 

3) (a handle of index 2) to a tubular (ribbon) neighborhood of a circle 
in the boundary of N^M°'~^ it is glued a thickened two-dimensional disc 
therewith the boundary of disc is glued to the circle. 

We call a critical point degenerate otherwise. The number i is called the 
index of a critical point. 

The following theorem is proved straightforwardly by the enumerating 
combinatorial types of critical points. 

Theorem 1 All possible type of critical points are given in Tahle\^ 

The types 1-14 correspond to nondegenerate critical points. The last type 
correspond to the "monkey saddle" which in the smooth case if given by a 
critical point (x, y) = (0, 0) of the function f{x, y) = x'^ — xy"^ . The passing 
through this value results in gluing two handles of index 1 . 

Let Mj, i = 0, 1, 2, be the sets of critical points of index i. Every monkey 
saddle v contributes to Mi the two points: the point v itself and its formal 
"double" v' . For every i = 0, 1, 2 let us consider the vector space Cj over Z2 
with generators from Mj. To construct a chain complex we have to define 
the linear boundary operators 9i : Ci — )• Co and ^2 : C2 — > Ci. 
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Table 1: Types of critical points 
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3 The gradient flow and the boundary operators 



The gradient flow enables us to describe a gluing of handles and that allows 
us to compute the boundary operators in the chain complex C*. By a 
trajectory of the gradient flow we mean a sequence of vertices of M which 
we may pass going down along edges. The difference from the gradient flow 
of a smooth function consists in an indeterminacy, i.e., from every vertex 
we may descend into not a unique vertex that lies below. Since we use 
the Alexander duality, it is enough to define the gradient flow from critical 
points of index 1 to critical points of index 0. 

Let us recall that we compute homology with coefficients in the filed Z2 
which consist of the two elements and 1 meeting the summation rules: 
+ = 1 + 1 = 0, + 1 = 1 + = 1. 

By a tangent vector ^ = {v,e) € TM we mean a pair consisting of a 
vertex v and an edge e 3 v which lie in the boundary of M. We say that 
^ is increasing if another end of e lies at a higher level of /; otherwise, we 
say that C is a decreasing vector. Moreover on the set of increasing (or 
decreasing) tangent vectors at a given point we fix a natural order that is 
induced by ordering of the coordinate directions. 

Let w be a critical point of /. If v has index 0, then exactly three 
tangent vectors are drawn from v and all of them are increasing. If u is a 
nondegencratc point of index 1, then to every edge of a decreasing tangent 
vector we assign some label S that is equal to ±1 and is as follows: The list 
of all types of critical points shows that all such edges lie in the boundary 
of M and there are at least two of them. If there are exactly two such edges 
then we mark one by +1 and another by —1. If there arc three such edges 
then two of them belong to an elementary face from the boundary of M and 
we mark these two edges by the same label, say by +1, and then mark the 
left edge by —1. 

Monkey saddle. If f is a monkey saddle, then wc add a fictive 
critical point v' , which lies above v with respect to the level surface of 
/, and also add a fictive edge vv' , connecting v and v' . We put that 
{v', (1, 0, 0)), («', (0, 1, 0)), {v, (0, 0, 1)) and {v, vv') are increasing vectors and 
(f', (0, 0, — 1)), {v',v'v), (f', (— 1, 0, 0)), and (0, — 1, 0)) arc decreasing vec- 
tors. Wc put that the end of every non-zero vector (v' , (a, f5,'~f)) coincides 
with the end of {v, {a, /3, 7)), a, /3, 7 = —1, 0, 1. 

Let t; be a critical point of index 0, and v' be a critical point of index 1. 
A pair of tangent vectors ^ = {v, e) and ^' = {v', e') is called a gradient 
pair, if ^ is an increasing vector, ^' is a decreasing vector, and the following 
conditions hold: there exists a sequence of increasing tangent vectors (e^, fj). 
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i = 0,...,k such that cq = e, vq = v, the edge Cj connects Vi and v^+i, 
Cjt = e', and for every vertex Vi the vector {ei,Vi) is greater (with respect to 
the order) than ah other tangent vectors which start at this vertex and he 
in M. 

Exactly in this case we will write that ^ < We mean by the gradient 
flow of / on M the family of gradient pairs of tangent vectors: 

GF{M) = m,^2m eTM,i = 1,2,6 < 6}- 

By definition, a path in M is a family of successively traversed edges from 
M. The path formed by edges eo, • • • , from the definition of a gradient 
pair is called gradient. 

We define the boundary linear operator by the formula: 

divi = v^ + Vq, Vie Ml, 

where 

v+ e{ve Mo|3e = iv,e),3^i = ivi,ei),S{ei) = +1,^ < 6}, 
voe{ve Mo|3^ = {v,e),3^i = {vi,ei),S{ei) = -1,^ < 6}- 

We recall that Mq and Mi are the sets of critical points of indices and 1. 

Let us remark that this definition is correct because from every critical 
point of index 1 there are drawn two decreasing tangent vectors labeled by 
different signs. Therewith we may go down until we achieve (in finitely 
many steps) a critical point of index 0. Moreover Vq and Vq are uniquely 
defined because we always choose the highest tangent vector. If we consider 
two gradient paths 71 and 7^ going down to Vq and Vq , then together these 
paths form a path whose neighborhood forms a handle of index 1 attached 
when the value of / passes through f{vi) (see §2). 

Let us consider the three-dimensional body M' = K\M, the complement 
to M, and a function h = -f on M' . Clearly, M'- = M2-i, i = 1,2; 
Mq = M2 U {po}, where po is the corner point of K, where h achieves its 
minimal value. Let C^' be the vector space over Z2 generated by the elements 
of M-. We have the isomorphisms C'^ = C2~i, i = 1,2, and Cg = C2 Z2, 
where the Z2 component is generated by pQ. By applying the construction 
of di to h we obtain the boundary operator dih : C[ ^ C'q. Let us consider 
the dual spaces C'*, i = 0, 1,2, consisting of Z2-valued linear functions on 
C[. There are isomorphisms C* = C[, i = 0,1,2, which assign to every 
Z2-function on C[ the formal sum of points from M[, at which the function 
does not vanish. The operator dih defines the adjoint operator 

5q = d*ih : = C2 ® Z2 ^ C'l* = Ci. 
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We put 

^2 = 5olc2 • ^2 — ^ C*!- 

Theorem 2 The homology groups of the chain complex C* = {Ci,di) are 
isomorphic to the homology groups of M with coefficients in Z2 : Hi{Cif) = 
H,{M;Z2), i = 0,1,2. 

Proof For every critical point of index 1 there exist two paths such 
that they go down to critical points of index in the directions labeled by 
different marks. This pair of gradient paths is combined into the path that 
connects two critical points of index and hence the formal sum of these 
points belongs to the image of di . By using the standard reasonings of Morse 
theory, it is easily proved that if two critical points of index lie in the same 
connected component of M then they are connected by the chain of such 
pairs of gradient paths. Since every connected component of M contains at 
least one critical point of index 0, we have Hq{M) = Co/lmdi. Analogously 
it is shown that H^{M') = Ker^o- Then we have H^{M') = Keidp,. Here 
we recall that the reduced homology groups Hi meet the equalities Hi{M) = 
Hi{M) for i > and Hq{M) = Hq{M) © Z2; the analogous equalities are 
valid for the reduced cohomology groups. 

Let us recall the Alexander duality theorem (see, for instance, |13j): 
Let A C M" be a compact polyhedron. Then we have the isomorphism: 

F,(M"\A;Z2) «F"-5-i(^;Z2). 

By definitions and duality, we have H2{M) = H2{M) = H°{M') = 
Ker {82). Hence the 0-th and the 2-nd homology groups of C* are isomorphic 
to the corresponding homology groups of M. 

The Euler characteristic xi^) of M is equal to the alternated sum of the 
numbers of cells of any cellular decomposition of a deformation retract of M 
and, in particular, to x{C^), and it is also equal to the alternated sum of the 
Betti numbers, i.e. to 60 — + &2- Since x{C*) = &o(C'*) — hiC*) + 62 (C**), 
we have 

X(M) = x{C,). 

It suffices to notice that 61 (M) = 6o(M) + 62(M)-x(M) = 6o(C*)+fe2(C*)- 
x{C*) = 61 (C*), and, since we consider homology with coefficients in a field, 
Hi{M) = ifi(C*). The theorem is proved. 
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4 An algorithm for computing the boundary ope- 
rators 

The approach described above may be used for constructing an algorithm as 
follows: Initially we have some array M = M[i, j, k] , where i,j,k = 1, . . . ,N, 
whose entries are equal to or 1. Therewith the three-dimensional body 
M is formed by all elementary cubes + 1] x + 1] x [k,k + 1] with 
M[iJ,k] = l. 

The preprocessing consists in the following. We stretch the array M 
in every direction, i.e. replace it by the array M' such that M[i,j,k] = 
M'[3i + h,3j + ji, 3k + ki], ki = 0, 1, 2, i,j,k = 1, . . . N. Then we 

assign zero values to such elements M'[i,j, k] that correspond to elementary 
cubes [i,i + 1] x + 1] x [k, k + 1] intersecting the boundary of M' over 
a two-dimensional face. 

At the input of the algorithm we have a preprocessed array M. At the 
output we obtain the lists CO and CI of critical points of indices and 1 (as 
it is described above every monkey saddle contributes into CI two elements; 
these lists give bases for the groups Co and Ci) and the matrix D = D[i,j] 
that defines the boundary operator di. 

To every vertex v we correspond a formally adjoint vertex v' . Let us 
consider the function Ind which assigns to every vertex v € M either its 
index, if the vertex is a critical point of index or 1; either —1, if it is a 
monkey saddle; either —2 otherwise. The list CO contains the critical points 
of index 0, and the list CI consists of the critical points of index 1 and the 
points V and v' which correspond to every monkey saddle v. The function 
Num returns the number of a critical point of index i in Ci. To every 
vertex v there corresponds GF{v), the number of a certain critical point 
of index to which one may go down from v. Finally, to every vertex v 
the functions DownA and DownJl return the ends of decreasing tangent 
vectors starting at f , where DownA corresponds to the highest vector, or the 
empty values if there are no such vectors. Therewith to the critical points 
of index 1 the functions DownA and DownJl return vertices corresponding 
to decreasing vectors with different signs. If f is a monkey saddle, the we 
assume that DownA{v) and DownJl{v) are the ends of vectors (—1,0,0) 
and (0, —1, 0) attached to v. For the adjoint vertex we put DownA{v') = v 
and Down_2(v') is the end of vector (0, 0, —1) attached to v. The description 
of the algorithm is given in Table [2j 

Applying this algorithm twice to M and its complement M' , obtain the 
matrices Dl and D2 describing the operators di and 82- Afterwards the 



10 



Table 2: Algorithm 
Input: the preprocessed array M 

Output: the lists of critical points CO and CI, the boundary operator 
matrix D 



1: D := 0, CO = 0, CI = 

2: for for all integer- valued vertices v e M do 

3: if Ind(v) = then 

4: Add{CO,v) 

5: GF{v) := {Num{v)} 

6: if Ind{v) = 1 then 

7: Add{Cl,v) 

8: GF{v) := GF{Down_l{v)) 

9: if GF{Down_l{v)) / GF(£»ou'n_2(?;)) then 

10: D{Num{v),Num{GF{DownA{v)))) := 1 

11: L>(A^«m(u),iV'um(GF(Dou;n_2('y)))) := 1 

12: if Ind{v) = -1 then 

13: Add{Cl,v) 

14: Add{Cl, v') 

15: GF(7;) := GF{DownA{v)) 

16: GF(i/) := CF(Doii;n_l(f)) 

17: if GF(Z)ou;n_l(?;)) / GF(£»ow;n_2(?;)) then 

18: D{Num{v),Num{GF{DownA{v)))) := 1 

19: D{Num{v),Num{GF{Down.2{v)))) := 1 

20: if GF{DownA(v')) / GF(Dou;n_2(t;')) then 

21: D{Num{v'),Num{GF{DownA{v')))) := 1 

22: D{Num{v'),Num{GF{Down.2{v')))) := 1 

23: if Ind(v) = ~2 then 

24: GF{v) := GF{DownA{v)) 

25: end 
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Betti numbers are found from the formulas: 

bo = dim (CO) - rank (Dl), 

bi = dim (CI) - rank (Dl) - rank {D2), 

b-2 = dim (C2) - rank (i:i2). 

The complexity of the algorithm for constructing the matrices of the bound- 
ary operators, including the preprocessing of data, is equal to 0(n) where 
n = is the number of vertices in K. For a subsequent computation of 
Betti numbers we have to compute the ranks of Di, for instance, by the 
Gauss method, and that needs 0(n|!) operations where is the number of 
critical points. Generically the relation between n and ric depends on the 
initial data. 
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